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Abstract. The performance of the positive P phase-space representation for exact 
many-body quantum dynamics is investigated. Gases of interacting bosons are 
considered, where the full quantum equations to simulate are of a Gross-Pitaevskii 
form with added Gaussian noise. This method gives tractable simulations of many- 
body systems because the number of variables scales linearly with the spatial lattice 
size. An expression for the useful simulation time is obtained, and checked in numerical 
simulations. The dynamics of first-, second- and third-order spatial correlations are 
calculated for a uniform interacting ID Bose gas subjected to a change in scattering 
length. Propagation of correlations is seen. A comparison is made to other recent 
methods. The positive P method is particularly well suited to open systems as no 
conservation laws are hard-wired into the calculation. It also differs from most other 
recent approaches in that there is no truncation of any kind. 
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1. Introduction 

Simulation of exact many-body quantum dynamics is a highly non-trivial task even in 
the simplest models. However, such a capability is highly desirable for studies of the 
quantum behaviour of mesoscopic systems. First principles methods are particularly 
useful when there are several length, time, or energy scales of similar size, which makes 
simplifying approximations inaccurate. It is well known, however, that a calculation 
using an orthogonal basis is intractable because the size of the Hilbert space needed 
to accurately describe the state grows exponentially with the number of subsystems. 
Perturbation theory is often not useful either, owing to questions of convergence. Other 
commonly used techniques - like the mean-field approximation - rely on factorizations 
of operator products, with unknown regimes of validity. 

| www.physics.uq.edu.au/BEC 
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A promising approach is to simulate stochastic equations derived from phase-space 
distribution methods such as the positive P[TJ |2] and gauge P[3| representations, or 
representations based on other separable basis sets. In such methods the quantum state 
is written as a probability distribution of off-diagonal system configurations, which are 
separable between local subsystems. For the case of the positive P representation, these 
subsystems are the spatial lattice points. The number of variables needed to specify a 
single configuration grows only linearly with the system size because they are separable. 
Correlations between subsystems are contained in the details of the distribution, which 
is stochastically sampled. This very mild growth of the number of required variables is 
the reason that such methods lead to tractable many-body simulations. Recent work 
has shown that any model with two-body interactions can be expressed by such a 
distribution of separable system configurations^]. 

In this paper we develop reliable ways of estimating the feasible simulation time, 
which is limited in practise by the growth of the statistical error. This allows an 
assessment of simulation parameters before a potentially lengthy calculation is started. 
We apply our method to single and coupled anharmonic oscillators, and compare 
analytic and numerical results for the simulation time. We also investigate the quantum 
dynamics of a one-dimensional interacting Bose gas. This calculation shows a unique and 
interesting behaviour under conditions of a sudden switch in the interaction strength. 
We observe correlation waves in which an otherwise homogeneous quantum gas develops 
wave-like structures only found in the relative coordinates. 

This is a particularly timely issue in view of recent developments in fibre optics 
and ultra-cold atomic Bose gases. It is now possible to experimentally investigate 
large numbers of interacting bosons, under conditions where quantum coherence plays 
an important role. As an example, quadrature squeezing in quantum solitons was 
predicted using these techniques, and observed experimentally. Ultra-cold atomic Bose 
gas dynamics has been observed in a number of geometries. In these experiments, it is 
possible to modify external potentials in order to drive the quantum gas into a state 
very far from either the ground state or any thermal equilibrium state. More recently, 
experiments have demonstrated the possibility of dynamically varying the inter-particle 
potential. 

The main body of the paper is organized as follows: Firstly the positive P method 
is summarized in Section and briefly compared to other approaches in Section El 
Subsequently in Sections Eland EH the estimates of useful integration times are obtained, 
and their implications considered. These estimates are verified by numerical simulations 
of single-mode and multi-mode systems, in Sections El and respectively. As an 
example, in the latter Section, the dynamics of spatial correlation functions in a ID 
uniform gas are investigated. We consider the effects of an experiment in which the 
inter-particle potential is suddenly changed. This is potentially observable using the 
technique of a Feshbach resonance in an external magentic field. 

In a subsequent paper we consider the more sophisticated stochastic gauge 
method for these types of simulations. 
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2. Model 



2.1. The lattice model 

We consider a dilute interacting Bose gas. In the continuum, the second- quantized 
Hamiltonian with kinetic terms, two-body interactions with interparticle potential U(r), 
and external potential Kxt( x ) is 



H = 



^V$t( x ) V $(x) + KxtW^x^x) 

+ Pr ( X - y ) $t( X )$t(y)$( X )$(y) 



(1) 



d 3 x. 



Here ^(x) is a boson field operator at x and m is the boson mass. 

In practice, (JTJ is replaced by a lattice Hamiltonian, which contains all the essential 
features provided the lattice spacing is sufficiently small. For a rarefied gas of the kind 
occurring in contemporary BEC experiments, s-wave scattering dominates jH], and the 
s-wave scattering length a s is much smaller than all other relevant length scales. If the 
lattice spacing is also much larger than a s , then the two-body scattering is well described 
by an interaction local at each lattice point. 

Let us label the spatial dimensions by d — 1, . . . , D, and label lattice points by the 
vectors n = (m, . . . , n D ). For lattice spacings Ax d , the spatial coordinates of the lattice 
points are x n = (mAxi, . . . ,nz>Axz>). The volume per lattice point is AV = Yl d Ax^. 
We also define the lattice annihilation operators a n ~ vAK $(x n ), which obey the usual 
boson commutation relations of [a n ,aj 



4] 



<5nm- With these definitions, one obtains: 



H = h 



E 



1 



2 ^ 

n 



(2) 



In this normally ordered Hamiltonian, the frequency terms a; nm = oo^ nn come from the 
kinetic energy and external potential. They produce a local particle-number-dependent 
energy, and linear coupling to other sites, the latter arising only from the kinetic 
processes. The nonlinearity due to local particle-particle interactions is of strength 

with the standard coupling valuejH] being g = g 3D = 47ia s h 2 /m in 3D, and g = gz^jo in 
2D and ID, where a is the effective thickness or cross-section of the collapsed dimensions. 

When interaction with the environment is Markovian (i.e. no feedback), the 
evolution of the density matrix p can be written as a master equation in Linblad form 



dp 
di 



1 

ih 



H,p 



(4) 



For example, single-particle losses at rate 7 n (at x n ) to a T 
by L n = a^J^u- 



heat bath are described 
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2.2. Positive P representation 

The positive P representation was developed in |2j. Here we summarize the issues 
relevant to the dynamics of the model (JTJ), and present the stochastic equations to 
simulate. 

For a lattice with M points, the density matrix is expanded as 

p = J P+(a, (3) A(a, (3) d 2M cx d 2M (3, (5) 

where A is an off-diagonal operator kernel, separable between the M lattice point 
subsystems. We use a coherent- state basis so that 

A(a,/3) = || a){(3* || exp[-a-/3], (6) 

in terms of Bargmann coherent states with complex amplitudes a = (. . . , a n , . . .): 

1 1 a) = ® n exp [a n al] | 0). (7) 

The gauge P representation mentioned previously [3J includes an extra global weight in 
the kernel, which allows useful modifications of the final stochastic equations |3J 13 IE]- 
It has been shown that all density matrices can be represented by a positive 
real -P+J2J- The expansion (jSJ) then becomes a probability distribution of the A, or 
equivalently of the variables ex., (3. A constructive expression for P+(p) is given by 
expression (3.7) in reference [2], although more compact distributions may exist. In 
particular, a coherent state \a ) will simply have 

P + = 6 2M (a- ao )5 2M (f3- a :). (8) 

Using the identities 



a n A = a n A, (9a) 
() 1 A,, (96) 



LA 



/9„ + 



the master equation (j3J) in p can be shown to be equivalent to a Fokker Planck equation 
in P + , and then to stochastic equations in the a n and /3 n . The standard method is 
described in reference jS]. The correspondence is in the sense that appropriate stochastic 
averages of a n and j3 n correspond to quantum expectation values in the limit when the 
number of trajectories S — > oo. In particular one finds thatjZj 

n^>*) = (n +n^<) • (io) 

jk I \ jk jk I s 

Any observable can be written as a linear combination of terms (|10p. We will use 
the notation (-) s to distinguish stochastic averages of random variables from quantum 
expectations (•). 

The resulting 2M Ito complex stochastic equations are found using methods 
described in Refs. (21 E]. For the model (J2J) obeying the master equation (|3j) with 
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coupling to a zero temperature heat bath, they are of the Gross-Pitaevskii (GP) form, 
plus noise: 

da n =[-i ^2 w nmOm - inn n a n - 7 n a n /2]c/t + ^ B^} dW k , 

k (n) 

dPn=[ iJ2 ^nm/^m + ^n n /3 n - i n n /2]dt + £ B$> dW k . 

m k 

Here, n n = a n /3 n , and there are M' > 2M labels k to sum over. The dW k are 
independent Wiener increments (dWj(t) dWk(s)) s = 5jkS(t — s)dt 2 . In practice, these 
can be realized at each time step At by independent real Gaussian noises of mean zero 
and variance At. The elements of the M x M' noise matrices B must satisfy 

£ B iM= iK ' 5 »' 5 ™' (12) 

These do not actually specify the elements of B completely. The usefulness of these 
degrees of freedom are considered in detail in Refs. [SI El IE], but here we will just use 
the simplest choices: 

B ^k =iViKat nj 5 jtk , 

B n-k =V / ^/?n/(i+Af),fc 

(with k = 1, . . . , 2M). The indices j = 1, . . . , M label the lattice points. 
The simulation strategy is (briefly): 

(i) Sample a trajectory according to the known initial condition P+(0) = P + (p(0) ) 

(ii) Evolve according to the stochastic equations calculating moments of interest, 
and recording. 

(iii) Repeat for S 3> 1 independent trajectories and average. 
3. Comparison with other dynamical methods 

Several positive P and gauge P simulations of the models of Section ETTI (or very similar) 
have already appeared: 

• Quantum optical soliton propagation in nonlinear Kerr media[10j, which obeys a 
Hamiltonian formally very similar to 

• Dynamics of evaporative cooling and incipient condensation of an interacting Bose 
gasjmCE]. 

• Thermal and dynamical behaviour of spatial correlations in one-dimensional 
interacting Bose gases [TB^ 15] using the stochastic gauge represent at ion [7J Ej. 
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Some of the other many-mode dynamics simulation methods proposed more recently 
include (See Tabled for a visual comparison): 

The stochastic wavefunction method of Carusotto, Castin and Dahbard|14[ ITH] 
shares many common features with the gauge P method jH]- It is also based on a phase- 
space distribution over separable operators with a global weight, and can be regarded 
formally as a particular choice of gauge and basis set. In this case the separable 
"subsystems" are chosen to be exactly N identical particles, rather than the M spatial 
lattice points of the gauge P/positive P approach. Hence, sampled variables specify a 
orbital occupied by these N particles, rather than the conditions at each lattice point. 
This explicitly imposes both a precise particle number and particle conservation on 
the model, and so makes the method directly applicable only to closed models where 
processes such as random losses, outcoupling, or gain are absent. Where both stochastic 
gauge and stochastic wavefunction methods can be applied, the useful simulation time 
is similar [5J. 

A different approach entirely is based on extensions of the density matrix 
renormalization group (DMRG) met hod [151 ITT] . This method is based on the concept 
that an iV-subsystem state is equivalent to a construction called a matrix product 
state (MPS) in a larger Hilbert space. This involves a further N "ancillary" systems 
entangled with the original subsystems. In cases where entanglement between the 
original subsystems is small or short-range (such as spin chains with nearest-neighbour 
interactions), a good approximation to the full state can be obtained by truncating the 
amount of entanglement present in the MPS construction. These truncated states can 
be stored efficiently with a number of variables polynomial in N, due to the way the 
MPS states are constructed. This can lead to tractable deterministic calculations even 
with significant N. For longer-range coupling in 2D and 3D, there are indications based 
on the entropy of bipartite entanglement that the amount of entanglement is expected 
to grow too fast for tractable accurate simulations with this method[TH]. 

While each approach has its own merits, for large and/or 2D, 3D systems the 
stochastic phase-space methods appear to be the most efficient, perhaps apart from 
special cases. There, for the price of limited precision, one can obtain an un-truncated 
first-principles simulation whose demands on resources and time scale only linearly with 
size. For open systems, the methods based on lattice point subsystems (such as the 
positive P) can be applied in a straightforward manner, as they explicitly allow a varying 
particle number. In contrast, attempts to use orbital-based states as in the stochastic 
wavefunction lead at best to additional complication. 

4. Single-mode simulation times 

4-1. Exactly solvable single-mode model 

The single-mode case is of special interest as it is exactly solvable, while still describing 
an interacting many-body quantum system. It already contains the essential features 
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Table 1. Some general features of several methods for simulation of many-body 
dynamics. M is the number of lattice points, N the number of particles, d is the 
Hilbert space dimension of the subsystem at each lattice point, D the dimension of 
each corresponding ancillary state after possible truncation of entanglement in the 
DMRG-based methods. 

stochastic gauge stochastic DMRG-based 
positive P wavefunction methods 



subsystems 
stochastic 

hard-wired conservation laws 

open systems 

truncation 

tractable coupling 

memory requirements 

Hilbert space dimension 

references 



lattice points 
yes 
none 
yes 
none 
wide variety 
oc M 



oc 



M 



[U 121 13 El EH 13 



particles 

yes 

particle number 
no 
none 
wide variety 
oc M 
M N 



d— level systems 
no 

(i-level systems 
yes 

entanglement 
short-range or ID 
oc MdD 
d M 

[iniEEZlEEE] 



that make the distribution of trajectories broaden, thus limiting the useful simulation 
time. 

To simplify the notation, the mode labels n = (0, . . . , 0) will be omitted when 
referring to the single-mode system. Furthermore, we move to an interaction picture 
where the harmonic oscillator evolution due to the ftwcfia term in the Hamiltonian 
is implicitly contained within the Heisenberg evolution of the operators. Then, this 
"anharmonic oscillator" model has 

H=^a j < 2 a 2 , (14) 

and the equations to simulate are 

da = al—inndt — 7 dt/2 + iViftdWi], 

1— ( 15 ) 
dp = P[ iKndt-jdt/2 + ViKdW 2 ], 

with n = a/3. Note that by (JTUJ), the mode occupation (cfla) is (Re [n]) s . 

These equations can be formally solved using the rules of Ito calculus as 

logn(t) = logn(0) -jt + Viiifait) +iCi(*)L ( 16 «) 
loga(t) = loga(0) + (in — j)t/2 + i\fin( i i(t) — in f n(s)ds, (16b) 

where 







rW(t) 

Q(t) = dW (16 c) 

Jo 

are Gaussian-distributed random variables of mean zero. However, they are not time- 
independent: 

(Cj(t)(k(s)) s = S jk min [t, s] . (17) 
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In what follows it will be convenient to consider the evolution of an off-diagonal coherent- 
state kernel 

A = A(a ,A)), (18) 

with complex "particle number" uq = aoA)- This is because for any general initial state, 
each sampled trajectory will start out as an Ao with some coherent amplitudes. 

With initial condition A , analytic expressions for observables can be readily 
obtained. In particular, we will calculate the first-order time-correlation function 

GW(0,t)=ft(a) = (19) 

which contains phase coherence information. Normalizing by Tr[Ao], 



(20) 



When the damping is negligible, n real, and the number of particles is n ^> 1, 
one sees that the initial phase oscillation period is 

*osc = -sin- 1 — « , 21 

and the phase coherence time, over which |G^^(0, t)\ decays is 

The first quantum revival occurs at 

^revival = ~- (^3) 
4-2. Causes of growth of sampling error 

From (|4.1|) . logn(t) is Gaussian distributed, while loga(t) is also Gaussian-like, at least 
at short times when the effect of the nonlinear term is negligible. The variables a and 
(3 = n/a appearing in observable calculations behave approximately as exponentials of 
Gaussian random variables, at least over short times. Hence, the properties of this kind 
of random quantity can tell us much about the behaviour of the simulation. 

If £ is a Gaussian random variable of mean zero, variance unity, let us define 

v = v e a t (24) 
with real positive a. The even m moments of £ are 

77? ' 

lim (£{t) m ) s = Try- — : — rr, (25) 

5^oo XSV ; /& 2 m /2( m /2)!' v ; 

while the odd m moments are zero. This can be used to obtain 

lim (v) B = v = v exp ( (26) 



and 



lim var 



=1 = e° 2 - 1. (27) 
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Figure 1. Number of samples S mm required to obtain a single significant digit of 
precision in v, the estimate of the mean of the random variable v — . 



I.e. the standard deviation of v is linear in a for a < 1, but grows rapidly once a > 1. 
With finite samples S, the fluctuations of (v) s around v are usually estimated to be 



S 

using the Central Limit Theorem (CLT). Hence, the number of samples needed to 
achieve a given relative precision is 

S oc var [v/v] . (29) 

This is shown in Figure^ In particular, we see that simulations with reasonable numbers 
of trajectories S < O (10 6 ) only give useful precision while 

a 2 < 10. (30) 

In actual calculations of moments, from (jlO)) we see that common observables are 
estimated by averages (/) s , where / is polynomial in a and/or n. Because the exact 
solutions (j4.1J) are Gaussian-like, / behaves similarly to the idealized variable v discussed 
above, and useful precision will be limited by 

var [log |/|]< O(10). (31) 

in analogy with (J3DJ). We have switched to O (10) rather than just 10 here, because the 
details of the distribution of log |/| may differ from a Gaussian at long times due to the 
nonlinear term oc J n in (|4.1j) . 

Observables that are first order in a, eft (such as G^'(0, t)), or first order in cfta will 
be limited by var [log |/| ] ~ var [log \ a(t)\ ) or var [log \n(t)\ ), respectively. 

4-3. Logarithmic variances 

The time evolution of the variances of the logarithmic variables can be calculated in the 
limit S —> oo from the formal solutions (|4.1|) . For the phase-dependent variables it is 
more convenient to consider just the mean of the variances of the two complementary 
variables a and j3: 

V = \ {var [ log \a(t) \ ] + var [ log \(3{t) | ] } . (32) 
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Consider coherent-state starting conditions p(0) = A /Tr[A ] such that n(0) = n for 
all trajectories. We will use the notation 

= n' + in'o, (33) 

etc. for real and imaginary parts. 

Noting that 1): logn is exactly Gaussian distributed, and 2): in the Ito calculus 
noises are uncorrelated with any variables at the same or previous times, one obtains 
from (jUTJ): 

V = Kt/2 + K 2 f ds [ ds'(n"{s)n"{s')) s . (34) 
Jo Jo 

Now if s' > s, then from ()16c|) (j(s') = Cj(s)+Q(s' — s), where Q(t) are independent 
of the Q but again with (Q(t)(k(s))s = <5jfcurin[t, s]. It is convenient to define 
C ± = [Ci ± (2]/V2. Then for s' > s 

(n"(s)n"(s% = e-^ +s '\e- 2 ^-^) s (e-^-^) s (cos[V<V - *)]>. (35) 

x {K) 2 (sin 2 [v^C + (^)])s + K) 2 (cos 2 [v^C + (^)])s + nX / (sin[2 v ^C + (^)])s} . 

The trigonometric averages can be evaluated using fl25(l . After some algebra and 
integration, one obtains 

1 -e-T*(l + 7t) 



lim V = Kt/ 2 — K 2 n' n >. 7 



1 2 1 |2 

+ k \n \ 



9-7 



1 - e~^ e~ qt - 1 
+ 



1 / 1 - e-T* 

2 V 7 



2 



(36 a) 



7 9 
where 

g = 2(7-«) (365) 
is a "damping strength" parameter. Also, 

var [log |n(i)|] = nt. (36c) 

4-4- Special cases 

Two cases are of particular interest: With no damping 7 = 0, V becomes 

V = Kt/2 - {Ktfn' /2 + -(Ktf\n \ 2 + O [{Kt) 4 \n \ 2 } (37) 
at low times 2nt <C 1, and 

V « Kt/2 + ^>o| 2 e 2Kt (38) 

at t > l/2«. 

If the damping is small but nonzero (^yt < 1), (|37|) becomes modified to 

V « «*/2 - n' (^) 2 [l - ~ 7 *]/2 + i|n | V) 3 [l - ^ 7 t] (39) 
to order (7/f:) 2 and (k£) 4 . Also, while q < 0, (J3S|) becomes V Kt/2 + /t 2 |n | 2 e~ IJ */g(g— 7). 
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With strong damping when q > 0, the long time (i.e. qt ^> 1) behaviour becomes 
asymptotic to a simple linear increase with time 

V = nt/2 - b, (40) 

where 

b=^[2n' -\n \ 2 e/(l-e)\ (41) 

is a constant, and e = k/j G [0,1]. We note b is positive for strong damping, but later 
tends towards — oo as e — > 1. The minimum damping required for positive b grows with 
mode occupation no. The short time behaviour when jt <C 1 is of the same form as the 
weakly-damped case above. 

Referring back to the condition (j31j) . we see that while the n(t) distribution remains 
usable for a relatively long time (t < O(10/k)) ~ l-5t rev ivah the phase-dependent 
variables can rapidly acquire very broad distributions (as per ()38jl) when damping is 
small and mode occupation is large. This is what limits simulation usefulness. 

4-5. Useful simulation time: analytic expressions 

While moments containing only the mode occupation n = a) a can be evaluated for long 
times of order t rev ivab this is due to special symmetries present only in the single-mode 
system. With inter-mode coupling, large spreads in any single a n will rapidly feed into 
the remainder of the variables via the coupling terms oc u; nm in the equations (flip. In 
particular, evolution of "occupation" variables n n is non-trivially coupled to the "phase" 
variables a n . With the M-mode systems in mind, then, we are primarily interested in 
precision of observable estimates in the single-mode model that involve phase variables 
a or P, not just n (e.g. G^ 1 ). 

Using the limit (J3TJ) and (|36ajl or its special cases in Section [Ol useful simulation 
times tgim can be estimated for coherent-state initial conditions 1)18)1 . For non-coherent- 
state initial distributions each separate trajectory still starts out as a sample of an 
off-diagonal coherent-state kernel ()18)1. The useful simulation time is then determined 
by the shortest useful simulation time to be expected from any of the S samples. 

For large mode occupation with no damping (the worst case), the term in 
\no\ 2 (nt) 3 dominates V at short times nt <C 1, and 

t- «-°®- (42) 

''Sim ~ I |9/Q ■ K^^J 

re|no| ' 

Checking back, nt ~ O (3)/|no| 2//3 , so ([42)1 is consistent with the short time assumption 
for | no | ^> O (5). At weak nonzero damping jt <C 1, (j4*2*)) becomes modified to 

W « O (3)/ K |n | 2/3 [l + O ( 7 / K |n | 2/3 )]. (43) 

At small occupations uq 1, simulation times are longer, and the regime 
where (l38)l applies is reached. For a significant range of |n | 1, the term oc |no| 2 
dominates V, and the useful simulation time scales very slowly with |no| as 

, , o(i)-|io g K| 



-q/A 



(44) 
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With weak damping, — q/2 — > k. At even smaller no, the nt term in V is greater, and 

i sim «O(10)/«. (45) 

For sufficiently strong damping such that q > 0, and gt S i m 3> 1, (f4f)J) applies, and 
long simulation times are possible: 

U « O (6 + 10)/«. (46) 
5. Empirical simulation time results 

To verify the analytic estimates, simulations of an undamped single-mode anharmonic 
oscillator were performed for a wide range of initial coherent states n = n' from 10 -5 
to 10 10 . The primary aim was to determine the actual dependence of useful simulation 
time t S i m on mode occupation, particularly for intermediate no values < O (1) where the 
simple expression (J42|) does not apply. This is shown in Figure 121 These results will be 
useful for subsequent many-mode simulation time estimates in Section El 

In what follows, the term useful precision for an observable O has been taken to 
indicate the situation where the estimate O = (f) s of (O) using S = 10 6 trajectories has 
a relative precision of at least 10% at the one sigma confidence level. This is assuming 
the CLT holds so that 



(47) 

is used to assess uncertainty in O. For the model here, we consider useful precision 
in the magnitude of phase- dependent correlations \G^(0, t)\, which is the low-order 
observable most sensitive to the numerical instabilities in the equations. 

Uncertainties in the calculated t sim times arise because the AIG^I were themselves 
estimated from finite ensembles of S = 10 4 trajectories. The range of t S i m indicated in 
Figure 121 was obtained from 10 independent runs with identical parameters. 

Taking the analytic scalings P^j) and (}44|) at high and low n into account, 
parameters in an approximate curve 

-1/C2 



t, 



est 








-1 


2 log 


( ^ + 0" 






\<« J. 





(48) 



have been fitted to the empirical data. Best values of Cj are given in Table El and 
in Figure 121 this fit is compared to the empirical data, and seen to be quite accurate. 
Ci characterizes the pre-factor for high n' , C3 a constant residual t smi at near vacuum, 
C4 characterizes the curvature at small n' , while ci is related to the stiffness of the 
transition. The expression (}4*%)) reduces to Cin' ~ 2//3 / k and (03 — C4 log n' )/ (k/2), when 
n ^> 1 and n' 1, respectively. Uncertainty Acj in parameters Cj was worked out 
by requiring E„ {[*est( c i ± A Ci ,n ) - t sim (n )] / At s[m } 2 = £ no {l + ([t es t( c i^o) - 
simvuyj/ ^ sim 



r sim 

(n )]/At sim ) 2 }. 
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Figure 2. Useful simulation time i s i m for a positive P simulation of undamped 
one-mode system 1)140 (solid lines) Triple lines indicate observed ranges. Dashed line 
indicates best estimate with 1)480 . 

Table 2. Empirical fitting parameters for maximum useful simulation time i g j m in 
one mode. The fit is to expression 1)480. 



Cl 


C2 


C3 


C4 


2.54 ± 0.16 


3.2 ± f 2 


-0.5 ± 0.3 


0.45 ± 0.07 



6. Multi-mode simulation times 

6.1. Extrapolation to many modes 

Firstly, let us note that for each mode n the nonlinearity that leads to rapid growth of 
distribution broadness at V ~ 10 depends only on the local variables a n and j3 n . This 
suggests that extrapolating the single-mode results of the previous section may lead to 
usable simulation time estimates for the coupled M-mode system. 

If the distribution of an amplitude variable in any mode acquires broad and/or 
rapidly growing tails (e.g. due to the single-mode instability caused by the nonlinearity), 
the coupling terms oc u nm rapidly disseminate its influence to the remaining modes, and 
all observables acquire large uncertainty. This effectively terminates the useful part of 
the simulation. From and/or Figure El we can see that for a single weakly-damped 
mode t sim drops off rapidly with increasing mode occupation n . At the same time the 
distributions of its phase variables broaden rapidly, and this broadness is distributed 
into the rest of the modes. Thus, the most highly occupied mode of a many-mode 
system limits the overall simulation time. 
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6.2. Simulation time estimate 

If we assume (for the time being) that inter-mode coupling does not significantly affect 
the time in which the most highly occupied mode reaches V ~ 10, one estimates that 

t sim ~ t est (max (t < tcstin) [n n (t)]), (49) 

where n n (t) is the mean occupation (a^a n ) at time t, and test is given by (|4*Kj) . Using the 
mean value n appears justified since for usefully precise simulations at higher occupations 
we expect \n(t) —n(t)\ <C n(t). For strongly-damped systems (j4*U|) applies rather than 

test- 

For a model arising as the discretization of a continuum field model (0) with density 
p(x), n n « p(x n )AU. When max(n n (t)) ^> 1, the estimated simulation time becomes 

t . „(^\ mm ( 5 o) 

What was perhaps not obvious before obtaining the above expression is that 
coarser lattices lead to longer simulations. The useful time scales as (AV) 1 ^ 3 
in the strong interaction limit. Thus, when simulating a field model, it is worthwhile to 
use the coarsest lattice for which physical predictions are not affected. 

6.3. Influence of kinetics in the multi-mode case 

Linear coupling between modes is what complicates the physical behaviour, and 
consequently also the noise behaviour in a simulation. Let us use the label z n for 
either a n or /3 n , and look at the short time fluctuations of the phase-space variables. 
These fluctuations can be separated as 

Szjfy « 5 dhcct z n (t) + 5 kinctic z n {t), (51) 

into the "direct" part due to local interparticle scattering as considered in detail for the 
single-mode anharmonic oscillator ()14j) . and secondary "kinetic" fluctuations. These are 
due mostly to transfer of direct fluctuations between modes by the kinetic inter-mode 
coupling, at least initially. 

If S direct z n (t) 3> # kmetlc z n (t), then properties based on the single-mode analysis (such 
as the simulation time estimates ()49|) and (jHUj) ) will be qualitatively accurate. A rough 
calculation for uniform gases at density p (see Appendix | Appendix A| for details) shows 
that at short times, 



var [\5 kinctic z n (t)\] _ 



h 2 t 2 ir 4 

(52) 



var [\5 dircct z n (t)\} [60m 2 (A\/) 4 / D _ 
while the right hand side is 1. This is independent of the uniform density itself, and 
indicates that the additional fluctuations due to kinetic effects: 

(i) Become relatively more important with time. 

(ii) Are more dominant for fine lattices. 
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Once (|52j) becomes > O (1), deviations from the simulation time estimates of Section lo^ 
can be expected, although this simple analysis does not tell us whether inter- 
mode coupling decreases or increases simulation time. In general a decrease seems 
intuitively more likely, but an increase may occur if d~ direct z n (t) and d~ kmetlc z n (t) become 
appropriately correlated to reduce overall fluctuations. 



7. Multi-mode example: ID gas 

7.1. Model 

As an example of many-mode dynamical simulations, we have simulated a uniform one- 
dimensional gas of bosons with density p and inter-particle s-wave scattering length 
a s . The lattice is chosen with a spacing Ax^ 3> a s so that interparticle interactions 
are effectively local at each lattice point, and the lattice Hamiltonian (j2J) applies. The 
stochastic equations to simulate are (JTTJ). Periodic boundary conditions are assumed. 

The initial state is taken to be a T = coherent wavefunction, which is a stationary 
state of the ideal gas (i.e. when a s — g — 0). Subsequent evolution is with constant 
a s > 0, so that there is a disturbance at t — when interparticle interactions are rapidly 
turned on. 

Physically, this would correspond to the disturbance created in a BEC by rapidly 
increasing the scattering length at t « by tuning the external magnetic field near a 
Feshbach resonance. Behaviour found in this uniform system will carry over without 
qualitative change to a ID BEC confined on a length significantly larger than the length 
scale of observed dynamical phenomena. Most of these phenomena extend over the 
order of a healing length 

e heal = -^=. (53) 

This is the minimum length scale over which a local density inhomogeneity in a Bose 
condensate wavefunction can be in balance with the quantum pressure due to kinetic 
effects [in]. 

A useful time scale here is 

/thealX2 h 

€ h 2pg K J 

(the "healing time"), which is approximately the time needed for the short-distance 
q ^heai-j i n ter-atomic correlations to equilibrate after the disturbance [8J. 



7.2. Spatial correlations 

Here we calculate the dynamics of spatial correlation functions, which has not been 
previously investigated from first principles in this model. We consider values averaged 
over the location of one particle out of the pair /triplet since this is a uniform gas: 
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First order: 

m y {a m a m ) {a m+n a m+n ) 

This is a phase-dependent correlation function. Its magnitude \g^'\ tells one the 
degree of first-order spatial coherence, while its phase gives the relative phase of the 
wavefunction at distances x n . 
Second order: 

— (2)/ x \ _ \ r ( a m Q m+n Q mQm+n) (56) 

M m (aLa m ) (aL+nflm+n) 
This describes the likelihood of finding two particles at a distance x n from each other, 
relative to what is expected of a coherent field. For a bunched field, g( 2 \x.) > 1 (e.g. a 
thermal state has ~g^(x) = 2), while spatial antibunching is evidenced by g^{x) < 1. 
Third order: 

—(3)/ \ _ 1 \ r (ftmQ m +n Q m-n Q m Q m+n Q m-n) fK7\ 
m \^m a m/ \ a m+n^m+n/ \^m— n^m— n/ 

This three-particle correlation function describes the likelihood of three particles with 
(equal) distances x n between nearest neighbours (relative to a coherent field). In a 
BEC, the rate of three-body recombination, which can limit the condensate's lifetime, 
is proportional to g^ 3 \0) 20J. 



7.3. Comparison to simulation time estimates 

A number of positive P simulations of a one-dimensional gas were made, with densities 
of p = 100/£ heal ,10/£ heal , and l/£ heal , and various lattice spacings. The number of 
lattice points was varied between simulations taking on values from M = 25 to M = 500 
depending on circumstances to encompass and resolve the phenomena seen. More lattice 
points were required at the higher densities where correlations were weaker, resulting in 
longer simulation times overall. 

Figure El compares the actual useful simulation times in simulations with what was 
predicted by expressions and (j5U|l . We can see that the fit is remarkably good 
despite inter-mode coupling. It is clearly seen also that the simulation time decreases 
with Ax = nj p. 

Note that values of t S i m are based on precision in g^ 2 \ Precision tends to become 
worse as higher-order moments of a or a) become needed, in e.g. g( 3 >. This is simply 
because higher-order moments require the averaged quantity / to be a higher-order 
polynomial in a and (3, and in effect var [log |/|] ~ c 2 V, with c being approximately the 
order of moment. 

Figures 0| to show calculated correlation functions for the example system 
p = l/£ heal . We see enhanced pairing/tripling of particles at several preferred relative 
distances. These preferred interparticle distances increase with time. Note that the 
local density p(x, t) is completely invariant, and this wave-like behaviour is seen only 
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n = pAx 



Figure 3. Comparison of useful simulation time i s j m to the estimates (|50|) (dotted) 
and (|49() (grey line). i s ; m based on precision in g( 2 \x). Data points for simulations 
of ID gases of densities 100/£ hcal , 10/£ hcal , and l/£ hcal , are shown with CIRCLES, 
triangles, and SQUARES, respectively. S = 10 4 trajectories. 




in the correlations. It is interesting to note that the overall disturbance advances at 
a much faster rate than the peaks/troughs of the correlations. We also see that while 
long-range coherence is steadily lost, there is a short distance scale corresponding to 
the distance of the nearest two-particle correlation peak on which phase coherence is 
enhanced compared to the long-range situation. 

Significantly more detail regarding calculations with the uniform gas and the 
correlation waves seen can be found in [Hj, Chapter 10. 
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Figure 5. First-order correlations. Details as in Figure^] 




Figure 6. Third-order correlations. Details as in Figure 01 



8. Conclusions 

We have obtained the following: 

• Simulation time estimates and (joTIj) for many-mode positive P simulations, 
the former using the fitting curve PKJ) and empirical data of Table El the latter 
applicable at maximum mode occupation ^> 1. It was found that useful simulation 
time decreases with density and increases with lattice spacing (coarseness). 

• These simulation time estimates were found to be accurate for simulations of 
uniform ID gases. 

• Dynamics of spatial correlation behaviour (1st to 3rd order) in an interacting 
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uniform ID gas after disturbance due to change in scattering length. Propagating 
"correlation waves" in 2nd- and 3rd-order correlations (Figures 0] to EJ) were seen. 

• Detailed characterization of the behaviour of a single spatial mode (jl4j) in an 
interacting Bose gas model. 

• Estimates ()52|) of relative importance of direct and coupling-induced fluctuations 
in phase variables at short times. 

A comparison with other more recently proposed methods was made in Table ^ and 
it is seen that each has its own merits, depending on the model. In broad terms 
the advantages of the positive P /gauge P method are its simplicity and versatility: 
Equations to simulate are just the GP mean-field equations with added Gaussian noise 
which, remarkably, is sufficient to account for all quantum mechanical features. Also, 
there is no truncation of entanglement, no hard- wired conservation laws, and open 
systems are the standard. The method is trivially adaptable to parallel computation: 
just run each trajectory separately. 

The dynamical timescales that can be simulated in a given system using the positive 
P method, can be easily assessed "on the back of an envelope" using the simulation time 
estimates given here. 

Appendix A. Scattering and kinetic fluctuations 

The aim here is to assess the relative strength of fluctuations in phase- dependent 
variables due to direct noise terms oc a/k and to mixing of noise (cx u) nm ) from other 
modes by kinetic processes. A uniform gas of density p, with 7 = V ext = is considered, 
and we use the notation of Section 16.31 We denote the linear coupling terms in the 



First-principles quantum dynamics in interacting Bose gases I 



20 



equations (jTTj) as 

dz n = - i w nm^m dt + . . . 
m 

= -i£^dt+..., (A.l) 

where z can denote either the variable a or (3*. Kinetic terms in the Hamiltonian ([T|l 



lead to 



— _ y |V_|2pik 5i -(x n -x m ) / a 9 \ 



on a lattice with M points, n = (ni, . . . , n_o) labels lattice points in a Fourier space 
such that kg = (Aki^i, • • • > Ak D N D ), where Ak d = 2n/L d , and L d is the full lattice 
extent in the d dimension. 

Since the gas is uniform, mean conditions at each lattice point must be identical, and 
random variables can be written as the sum of a global mean part and local fluctuations: 

£L z \t)=S (z \t) + 5S^(t) 
z a (t) = z{t) + 5z n {t) 



Due to lack of external forces \z(t) | = y/pAV. 

Let us denote the z n fluctuations when no inter-mode coupling is present (i.e. 
£n = 0) as 5 dircct z n . In this case the formal solutions (j4.1|) apply to each mode 
independently. At short times nt <C 1 the nonlinear drift term initially acts mainly 
as a deterministic time-varying phase, and so the dominant fluctuating contribution 
would be due directly to the noise terms. Then 

loga n (t) Psloga(i) + iViK(i n (t) 

log p u (t) « log p(t) + v^n it). 



with Q n being random variables of the (jl6c|) kind and independent for each n. At short 
enough times that var [log |z n (i)|] <C 1 , one then has 

z a {t) « z(t) + z(t) [log z a (t) - log z(t)} 

= z{t) + 5 diTCCt z n {t). (A.5) 

Using (|A.4|) . ()A.5|) . and the properties (JHJ), one finds 

var[|5 direct ^(t)|] (A.6) 

Now for the kinetic-induced fluctuations. At short times while these are small, i.e. 
var [J 5 kinctlc | ] <C var [| <5 direct | ] , the additional fluctuations in the local terms due to 
the mode-mixing can be ignored, and so the remaining fluctuations due to coupling are 
just 5 kmetlc z n (t) ~ —i Jq 5£n \s) ds. Using (|A.2|) . and substituting in for z m with the 
approximate (direct noise only) short time expression z m ~~z + <5 dircct ,2 m , one obtains 

X J2 (^ dirCC ^m'(^)^ dirCCt ^m(t / ))s|kfi| 2 |k^| 2 e -* a /-(x n -x m ,) e ik K .(x n -x m ) 



m,m',n,n' 
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Since the direct noise at each lattice point in the locally-interacting model is 
independent, then (5 dhect z* m ,{t")5 dhcct z m {t')) s = 5 mtm >pgmm[t" ,t']/h, similarly to (jPjl . 
After performing the integrations over t' and t" and simplifying the Fourier transforms 
one obtains 

hH 2 

Ylm 2 M 



var [|^ kinetic ^n(t)|] « var [\5 di ^z a (t)\] ^^^^|k fi | 4 . (A.8) 



For a .D- dimensional system with many modes, one can approximate 

MAV ! \\c\ 4 rl D h - C l( kmaX ) M?r4 

J2^ J |k| d 5{AV)±/ D 



I k *l - tt^jj / |k| 4 d-k = ^ AT ^ /n , (A.9) 



where ci(k max ) is a shape factor O (1) that depends on the ratios between momentum 
cutoffs fc™ ax = n/Axd in the various lattice dimensions. For example in ID, c\ = 1, 
while for 2D and 3D when the momentum cutoffs in each dimension are equal, one has 
ci = 3|, and c\ = 6|. Substituting into (|A.8|) one obtains the final estimate (JH2J) . 

References 



[1] 

[2] 
[3] 
[4] 
[5] 

[6] 



[7] 
[8] 



S. Chaturvedi, P. D. Drummond, and D. F. Walls, J. Phys. A 10, L187-192 (1977). 
P. D. Drummond and C. W. Gardiner, J. Phys. A 13, 2353 (1980). 
P. D. Drummond and P. Deuar, J. Opt. B-Quant. and Semiclass. Opt. 5, S281 (2003). 
M. Cetinbas and J. Wilkie, quant-ph/0407036 (2004). 



P. Deuar and P. D. Drummond, submitted to J. Phys. A: Math. Gen. (2004), jcond-mat/0501058 
See also HP- 

A treatment of these issues can be found in A. J. Leggett, Rev. Mod. Phys. 73, 307 (2001). 

For a more detailed treatment see J. Weiner, V. S. Bagnato, S. Zilio, and P. S. Julienne, 

Rev. Mod. Phys. 71, 1 (1999). 
P. Deuar and P. D. Drummond, Phys. Rev. A 66, 033812 (2002). 

In preparation. See also, P. Deuar, PhD thesis, The University of Queensland (2004), 



cond-mat/0507023 



[9] See e.g. C. W. Gardiner, Quantum Noise ( Springer- Ver lag, Berlin, Heidelberg, 1991); 

C. W. Gardiner Handbook of Stochastic Methods (Springer- Verlag, Berlin New York, 1983). 
[10] S. J. Carter, P. D. Drummond, M. D. Reid, and R. M. Shelby, Phys. Rev. Lett 58, 1841 (1987); 

P. D. Drummond, R. M. Shelby, S. R. Friberg, Y. Yamamoto, Nature 365, 307 (1993). 
[11] P. D. Drummond and J. F. Corney, Phys. Rev. A 60, R2661 (1999). 

[12] M. J. Steel, M. K. Olsen, L. I. Plimak, P. D. Drummond, S. M. Tan, M. J. Collctt, D. F. Walls, 

and R. Graham, Phys. Rev. A 58, 4824 (1998). 
[13] P. D. Drummond, P. Deuar, and K. V. Kheruntsyan, Phys. Rev. Lett. 92, 040405 (2004). 
[14] I. Carusotto, Y. Castin, and J. Dalibard, Phys. Rev. A 63, 023606 (2001). 
[15] Y. Castin, and I. Carusotto, J. Phys. B 34, 4589 (2001). 

[161 G. Vidal, Phys. Rev. Lett. 91, 147902 (2003); G. Vidal, Phys. Rev. Lett. 93, 040502 (2004); 

M. Zwolak and G. Vidal Phys. Rev. Lett. 93, 207205 (2004). 
[17] F. Verstraete, J. J. Garcia-Ripoll, and J. I. Cirac, Phys. Rev. Lett. 93, 207204 (2004); F. Verstraete 

and J. I. Cirac |cond-mat/ 0407066 (2004). 
[18] G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Phys. Rev. Lett 90, 227902 (2003). 
[19] See e.g. F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 71, 463 (1999), 
p. 481. 

20] Yu Kagan, B. V. Svistunov, and G. V. Shlyapnikov, JETP Lett. 48, 56 (1988). 



